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ABSTRACT 



< 



The use of the digital computer for the simultaneous generation of 
circuit responses and their corresponding sensitivity to parameter changes 
is considered for circuits containing diodes and transistors. Modeling of 
the diodes and transistors is based upon the Ebers-Moll equations, with 
the addition of voltage-dependent capacitors to account for switching 
time. The resulting general iterative equations are then used to calcu- 
late the effects of radiation on a p-n junction. The waveforms and 
circuit recovery time are studied as a function of circuit parameters. 

The use of sensitivity functions for predicting waveform variations is 
considered. 
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I. INTRODUCTION 



A measure of the sensitivity of circui t- response variables to param- 
eter variations is useful in circuit analysis and design. Many different 
definitions of sensitivity are available. H. Bode [1] defines a logarith- 
mi< 



ic or normalized sensitivity function of a variable, x.(t,a.)» with 

x. 9 In a. 1 J 

respect to a parameter, a . , as S^ 1 = - ?— . The inverse of this 

J a j d I M X^. 

definition is adopted by I. M. Horowitz [2]. Calculation of the sensi- 
tivity of linear-circuit responses to incremental changes in network 
parameters has been presented by J. Leeds [3]. Recently S. R. Parker [4] 

generalized the technique to include nonlinear circuits. Both Leeds and 

x. 9X. 

Parker define the sensitivity function as S = — — . This work was 



a j 



aa j 



preceded by that of P. Kokotovic [5] and R. Tomovic [6] who considered 
variations in the solution of a linear differential equation due to incre- 
mental changes in the coefficients. Tomovic defines the sensitivity 
x. 9x. . x i 

coefficients as S = — — , and Kokotovic adopts the definition S = 

S j “j 

8 1 n • , and discusses implementation of the sensitivity coefficient 

J 

using analog computation. The work by Leeds and Parker is oriented 

towards the digital computer. In this thesis digital calculation of 

sensitivity functions of nonlinear circuits is considered in detail with 

specific application to the study of the effects of radiation on solid- 

x. 9X. 

state junctions. The sensitivity function is defined as S = , = 

ax. a j 3 1n a j 

) • This definition is adopted since it enables sensitivity func- 

J J 

tions with respect to parameters having widely different nominal values to 
be conveniently compared on a percentage basis, and also avoids some 
difficulties encountered in the definitions adopted by Bode and Horowitz 
when the circuit-response variable is zero. It should be noted that the 
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above definitions of the sensitivity function apply only for infinites- 
imally small changes in x and a.. The practicality of this definition 

* J 

when applied to specific examples with large parameter increments is 
considered in this thesis. 

In Chapter II it is shown that for linear circuits the sensitivity 
function may be derived using the compensation theorem [7]. A convenient 
set of iterative equations is developed using the concept of state vari- 
ables which allows the time response of the states and their corresponding 
sensitivities to be obtained simultaneously. A linear example illustrates 
the procedure. 

In order to obtain useful results from the sophisticated circuit anal- 
ysis and design programs for nonlinear circuits, acceptable nonlinear 
models are required. Chapter III discusses modeling of diodes and tran- 
sistors for digital computer analysis based on the Ebers-Moll [8] equations. 
These models contain nonlinear capacitances and nonlinear resistances for 
the p-n junction A state-variable formulation is developed which results 
in a set of iterative equations whose solution yields time responses for 
circuits containing linear elements, diodes, and transistors. The proce- 
dure is illustrated by considering the effect of a radiation current pulse 
injected across (1) a p-n junction diode, and (2) the base-emitter p-n 
junction of a transistor operating in the active region. These examples 
yield new results concerning the effect of radiation on p-n junctions and 
comparative data for several runs is presented. 

Sensitivity analysis of -nonlinear circuits is discussed in Chapter IV. 
The state-variable formulation developed in Chapter III is extended, re- 
sulting in a set of iterative equations which allow time responses and 
their corresponding sensitivities to be obtained simultaneously. These 
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procedures are applied to the diode example of the previous chapter and 
the results interpreted in terms of the sensitivity of the diode response 
to percentage parameter variations. 
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II. SENSITIVITY ANALYSIS OF LINEAR CIRCUITS 



A. INTRODUCTION 

There are several general theorems which are useful in analysis and 
design of linear networks. One of these, the compensation theorem [7], 
has been known for many years and may be applied when it is desired to 
determine what effect an incremental impedance change in one branch of a 
network has upon the currents and voltages of other branches of the net- 
work. Recently the concept of the sensitivity of certain circui t-response 
variables to incremental changes in given network parameters has been 

A 

presented by J. Leeds [3]. This work was preceded by that of P. Kokotovic 

A 

[5] and R. Tomovic [6] who considered variations in the solution of a 
linear differential equation due to incremental changes in the coeffi- 

A A 

cients. Tomovic and Kokotovic implemented the sensitivity coefficient 
by using analog computation. Leeds' work is oriented towards the digital 
computer. 

In the next section sensitivity analysis of linear circuits is con- 
sidered. In the succeeding section the relationship between the older 
compensation theorem and the recent sensitivity concept is discussed. 

This correlation has not appeared in the literature heretofore. The final 
section contains an illustrative example to demonstrate the implementation 
of the results on a digital computer. 

B. DETERMINATION OF SENSITIVITY FUNCTIONS 
1 . Auxiliary Network Approach 

Leeds [3] developed a relationship between the sensitivity func- 

x i 3x i 

tion, S a = - ■ - ■ (where x. is a response variable and a. is a circuit 
a .i aa i 1 J 
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parameter), and node voltages and branch currents of an auxiliary network 
called the sensitivity model. The auxiliary network is obtained by re- 
ducing all independent sources to zero and adding a single source in the 
branch containing the element under consideration. The source depends 
upon the nature of the element. Fig. 2.1 summarizes the sources to be 
added for branches containing resistance, inductance, and capacitance 
respectively. The voltages and currents in the auxiliary network are the 
voltage and current sensitivities with respect to the parameter being 
considered. 

The advantage of Leeds' approach is that it allows the sensitivity of 
a network to be calculated by evaluating the responses of a coupled auxil- 
iary network rather than actually performing the differentiation indicated 
by the definition of the sensitivity function. 

2. Augmented State Equations 

Sensitivity functions may also be obtained using the concepts of 
state variables. The network is first represented in state-variable form 
as discussed by T. Bashkow [9], and Kuh and Rohrer [10]. 

x = Ax + Bu^ (2.1) 



where 

x is an n vector of circuit states 

A is a constant nxn matrix 

u^ is an m vector of forcing functions. 

x. ax. 

The sensitivity functions, S = — ^ , may be obtained by 

d • o in d* 

J J 

performing the indicated partial differentiation on (2.1). The resulting 
partial derivatives can be considered as additional states which may be 
solved in conjunction with the original state equations. They may be 
included in the following augmented matrix; 
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FIG. 2 . 1 . CIRCUIT PARAMETER BRANCHES AND 
SENSITIVITY MODEL BRANCHES. 
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If u_ i s independent of a,, that is if a. is a circuit parameter 

J %J 
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Separating the set of first-order partial differential equations 
of (2.3a) for each sensitivity parameter yields 

x. - Ax. + Bu + A . x + . u (2.5) 

-sj -SJ -S3 S3 - S 3 - v ' 

for j = 1 ,2, . . . , r 

Trapezoidal integration is commonly used to solve a set of first- 
order matrix differential equations because the procedure always converges 
and provides acceptable accuracy for most problems [11]. Utilizing trape- 
zoidal integration, (2.1) and (2.5) yield the following iterative solutions 

x(nT) = 4 >(T)x((n-l)T) + r(T) [u(nT)+u( (n-1 )T)J 

x sj (nT) = 4 >(T)x sj ((n-l)T) + r(T) [u, sj (nT)+u sj ((n-l )T)] (2.6) 

+ r 1j (T)[x(nT)+x((n-l)T)] + r 2j (T) [u(nT)+u( (n-1 )T)] 

for j=l ,2, ... ,r 
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where 



*(T) = [I - j-A]-'[I + | A] 

r(T) = [I - J-A ]" 1 ST 

T ,A,T 

r lj(T> * [I ' -4“ 

r y (T) = [I - jA ]' 1 



Equations (2.6) and (2.7) are a set of iterative equations whose 
solution yields the time response of the states and their corresponding 
sensitivities. Combining the results into a single equation results in a 
solution format similar to (2.3). 
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C. USE OF THE COMPENSATION THEOREM TO DERIVE SENSITIVITY FUNCTIONS 

The compensation theorem relates an incremental impedance change in 
one branch of a network with corresponding incremental current and voltage 
changes in the other branches [7]. In order to relate the compensation 
theorem directly with sensitivity analysis an analagous relation for incre- 
mental changes in the branch parameters R, L, and C is developed. 

1 . Resistive Change 

If the resistance in a branch is changed by an amount AR, and a 
voltage source equal to T(aR) is introduced (Fig. 2.2a) the network vari- 
ables remain unchanged. If another source is introduced (Fig. 2.2b) the 
branch voltage and current will change by an incremental amount a V and aI 
respectively. By applying the superposition theorem (Fig. 2.2c) the 
incremental voltage and current can be calculated. 

AV = AI • (R + AR) + I - CAR) (2.9) 

Dividing (2.9) by aR and letting AR-^0 results in (2,10). 

M - W R + 1 < 2 - 10 > 

2. Inductive Change 

If the inductance in a branch is changed by an amount aL, and a 
voltage source equal to aL • (^-) is introduced (Fig. 2.3a) the network 
variables remain unchanged. If another source is introduced (Fig. 2.3b) 
the branch voltage and current will change by an incremental amount a V and 
Al respectively. By applying the superposition theorem (Fig. 2,3c) the 
incremental current and voltage can be calculated, 

AV = (L + aL) ^ + aL ^ (2.11) 
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FIG. 2.3 INDUCTIVE CHANGES. 
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Dividing by aL and letting aL-> 0 results in (2.12). 



3 V _ i d / 3l \ , d I 

sT ' L dt dt 



( 2 . 12 ) 



3. Capacitive Change 

If the capacitance in a branch is changed by an amount aC, and a 
current source equal to aC-(^t) is introduced (Fig. 2.4a) the network 
variables remain unchanged. Introducing another current source (Fig. 2.4b) 
results in incremental changes in the branch voltage and current. Again 
by applying the superposition theorem (Fig. 2.4c) the incremental voltage 
and current can be calculated. 



Al = (C + AC) AC ^ 

Dividing (2.13) by aC and letting aC-> 0 results in (2.14). 



(2.13) 



3 1 = r d_ / 3_V\ , dV 
3C L dt 3C' dt 



(2.14) 



Equations (2.10), (2.12), and (2.14) are identical to the auxil- 
iary network equations of Fig. 2.1. This illustrates the direct 
relationship between the compensation theorem and sensitivity analysis. 



D. LINEAR EXAMPLE 

Consider the peaking circuit [12] in Fig. 2.5a operating in the linear 

region. The incremental equivalent circuit is represented in Fig. 2.5b 

where R] = r^ is the plate resistance and C is the coil capacitance plus 

output capacitance, stray capacitance and wiring capacitance. The equiv- 

^1 R 2 

alent circuit is indicated in Fig. 2.5 c where R = - 5 — r~5 ~ anc * v = 

R 1 R 2 R 1 H 2 

- r-tr^ yv i • 
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(b) 




(c) 



FIG. 2.5 PEAKING CIRCUIT. 
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It is of interest to determine how sensitive the output is with 
respect to incremental changes in R, L and C As representative values 
let L - 0.1 henry, C = 250yf arads , and R be adjustable according to the 
type of response desired. The response can be determined from the differ- 
ential equations of the circuit and is underdamped for K < 1 (R = 20), 
critically damped for K = 1 (R - 10), and overdamped for K > 1 (R = 5); 
where K = . 

The sensitivity functions may be obtained by constructing auxiliary 
networks as outlined in the preceding section, or by differentiating (2.1) 
with respect to In R, In L, and In C respectively. 

The state-variable representation is illustrated below. The solution 
takes the form of (2.8) where 

u^ = v = input voltage of 70 volts 
output voltage 
inductor current 
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The computer program used to implement the solution is included in 
Appendix B. This program is general in that it is easily modified to 
handle linear circuits containing more state variables and sensitivity 
parameters . 

The time response of the output of the circuit and the corresponding 
sensitivities are presented in Figs. 2.6 through 2.11 for underdamped, 
critically damped and overdamped responses. These results have been 
obtained using trapezoidal integration as well as a fourth-order Runge- 
Kutta integration technique. 
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FIG. 2.6. LINEAR EXAMPLE —UNDERDAMPED RESPONSE 
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SENSITIVITY FUNCTIONS (volts) 




FIG. 2.7 SENSITIVITY FUNCTIONS FOR UNDERDAMPED 
RESPONSE. 
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FIG. 2.8. LINEAR EXAMPLE" CRITICALLY DAMPED 
RESPONSE. 
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FIG. 2.9. SENSITIVITY FUNCTIONS FOR CRITICALLY 
DAMPED RESPONSE. 
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FIG. 2.10. LINEAR EXAMPLE -- OVERDAMPED RESPONSE. 
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FIG. 2.11. SENSITIVITY FUNCTIONS FOR OVERDAMPED 
RESPONSE. 
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III. MODELING OF NONLINEAR CIRCUITS 



A. INTRODUCTION 

The digital computer has played an increasingly important role in 
circuit analysis and design in the past ten years. This has been partly 
due to the emergence of sophisticated circuit-analysis and design programs 
[13]. In order for these programs to be useful, acceptable nonlinear 
models which relate model parameters to physical processes are required. 
These models must lend themselves to numerical analysis and provide accept- 
able quantitative accuracy. 

In the next section, a general mathematical format which is convenient 
for sensitivity analysis of nonlinear circuits is developed. In the suc- 
ceeding section modeling of diodes and transistors is outlined. This is 
followed by examples which consider calculation of the response of a 
p-n junction to a pulse of radiation current. 

B. GENERAL NONLINEAR STATE-VARIABLE FORMAT 

The state equation for a nonlinear time-invariant system is of the form 

x(t) = a(x(t) ,u_(t)) (3.1) 

where a is a vector of nonlinear functions of the state variables and 
forcing functions. 

As is shown in the section which follows, (3.1) takes the following 
general form for circuits containing diodes and transistors (represented 
by the Ebers-Moll [8,14] equations). This representation includes non- 
linear resistances and nonlinear capacitances. 

x(t) = A(x(t))*x(t) + c_( x ( t ) ) + B(xU))- u_(t) (3.2) 
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where A(x_(t)) is a nxn matrix and B(x_(t)) is an nxm matrix, both with 
elements which may be nonlinear functions of the states, and £(x_(t)) is 
a general nonlinear vector. 

Equation (3.2) may be written as 

x(t) = f(x(t)) + B(x_(t) )• u_(t) (3.3) 

where 

f(x(t)) = A(x(t))-x(t) + £(x(t)) 

Using trapezoidal integration for (3.3) results in 

x(n) - x(n-l) = ^ [f(x_(n)) + f(x(n-l))] 

+ Of [B(x/n))*£(n) + B(x(n-1 ) )-£(n-l )] (3.4) 

Equation (3.4) is a nonlinear matrix equation implicit in x_(n) which, 
as a result, is generally difficult to solve. Simplification results if 
f_U(n)) and B(£(n) )• u_(n) , where £(n) is considered constant, are expanded 
in a Taylor series about the point x_(n-l). In matrix notation 



f(x(n)) = z (ax^v)^ f(x(n-l ))/k! 

k=0 



B(x(n))»£(n) = z (ax^v)^ Cb (><_( n-1 ))-£(n)]/kl 

k=0 



where 



v = [ 



3 3 

3X-j 3X 2 



_3_] T 

3X n 



f, Bu, x, ax are n vectors 



and 



ax = x(n) - x(n-l ) 



(3.5a) 

(3.5b) 
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A first-order approximation for (3.5) results in (3.6). 

f(x(n)) = f(x(n-l ) ) + (ax T v) f(x(n-l)) (3.6a) 

B(x(n))- u_(n) = B(x(n-l))-u(n) + (ax T v) [B(x(n-1 ))*^(n)] (3.6b) 

Substituting (3.6) into (3.4) yields 

x(n) - x(n-l) - [ 2f(x(n-l)) + (ax T v) f(x(n-l)) ] 

+ ^-[B(x(n-1)) • (ujn) + u(n-l)) + (ax T v) {B(x(n-1 ))-u_(n)}] 

(3.7) 

Equation (3.7) may be simplified by using the matrix identity 

(ajJv) f = (vf T ) T ax (3.8) 

which is proven in Appendix A. 

Substituting (3.8) into (3.7) and solving for x(n) results in 
x(n) = _x(n-l) + aT- [I - ^-[(vf T ) T + { v[B”Ujn)] T } T ] ] _1 

[f + 1 b • (u(n) + u(n-l))] (3.9) 

where f, B, vf_^ and (B'U_(n)) T are evaluated at the point x(n-l). 

Equation (3.9) can be solved directly for xjn) from the known values of 
_x(n-l), u_(n), and u^( n-1 ) . 

> 

C. DIODE AND TRANSISTOR MODELING 

1 . Modeling for a p-n Junction Diode 

Diode modeling has received much attention during the past twenty 
years [15,16]. Recently Steele [17] prepared a report concerning semi- 
conductor modeling specifically for computer analysis. The applicable 
results are summarized in this section. 
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a. Current-Voltage Relationship 



The current-voltage relationship for an ideal diode is 




( 3 .. 10 ) 



where 



S 



reverse saturation current 



a = a parameter dependent upon the temperature and the 
typ.e of semiconductor material. 

Generally the reverse saturation current is not entirely 
independent of voltage. This is due to a leakage current and can be 
accounted for by the addition of a resistance, R^, shunting the diode. 

Also, at large reverse voltages an effect known as Zener breakdown occurs.. 
This effect is omitted in the following analysis. 

b. Diode Resistance 

The static resistance R^ ( v ) of a diode is defined as the ratio 

v/i^. Since R^v) varies widely with v and i^, it is generally not a 

useful parameter. An incremental, or dynamic, resistance r^v) is defined 
d v 

as r d (v) = ^p-— = — . Although the dynamic resistance varies 

d I,,ae a 

widely with v it is a useful parameter for small-signal operation. 

c. Capacitance Effect 

Diode capacitance is due to two separate effects; space charge : 
and charge storage. 

(1) Space Charge . A diffusion potential, V^, is established 
in the junction region when p and n materials are placed in contact. This' 
is caused by an increase in donor atoms on the n side and acceptor atoms 
on the p side of the junction, resulting in a capacitive effect. As a 
reverse voltage is applied more donor and acceptor atoms are uncovered 
causing the junction width to expand, and hence the capacitance to decrease. 
This voltage-dependent capacitance has the following form when the diode 
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is reverse biased. For radiation analysis, acceptable results are ob- 



tained when C D is assumed to be constant. 



C 



D 



(1 - v/v D ) n 



(3.11) 



where n = 1/2 for abrupt junctions 

n = 1/3 for li near-graded junctions. 

(2) Charge Storage . When the diode junction is forward biased 
the potential barrier is lowered. This results in a high density of minor- 
ity carriers injected into the p and n regions, which can be viewed as 
stored carriers. A change in voltage results in a change in density (num- 
ber of stored carriers) and hence a capacitive effect . 

If the voltage on the junction is abruptly reversed the 
stored carriers reach equilibrium after some storage time T^. The associ- 
ated current may be defined as 



i = T 



D 




(3.12) 



Substituting (3.10) into (3.12) yields 



T D i S ote 



av 



C (v) ^ 
L S l ; dt 



(3.13) 



where C<.(v) = I«-ae aV 1S the charge-storage capacitance. The charge- 
storage capacitance may be alternately written in terms of the dynamic 
resistance as C^(v) = T^/r^v) . 



d. Diode Model 

Combining the effects outlined above results in the diode 



model in Fig. 3.1. The circuit equation becomes 

v 



C d (v ^ dt + ] d + R 



0 



where C^v) - + C<-(v, 



SH 
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(3.14) 




Fig. 3.1. Diode model. 
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Equation (3,14) is in the form of (3,3) where 



f(x(t)) 






(3.15a) 



B ( x ( t ) ) = 0 



(3.15b) 



The diode example (III, D.) does not include parameter temper- 
ature dependence. The range of values studied includes a spread in 
parameter values as may be due to manufacturing tolerances, environmental 
conditions, and aging. The Zener breakdown effect and voltage dependence 
of Cq are neglected. Typical parameter values are used. 

2. Modeling for a p-n Junction Transistor 

Many models for a p-n junction transistor have been proposed in 
the last fifteen years, the most common being the Ebers-Moll model [8,14], 
the lumped model [16], and the charge-control model [18,19]. It has been 
shown that each model yields the same result in the solution of a large- 
signal transient problem [20]. The Ebers-Moll model is used in this study 
and is summarized in this section. 

a. Current-Voltage Relationship 

The current-vol tage equations for the Ebers-Moll model are 
given in (3.16) and (3.17). 



I 



E 



I 



C 



ot I I C0 

1 " 



(e aV C - 1) 



Veo 

1 - a i a N 



(e aV E - 1) 



EO 



a I a N 



L C0 



1 - aja N 



(e aV E - 1) 
(e aV C - i) 



where 

a = constant having a value between 1 and 2 
= common-base current in the normal mode 



(3.16) 

(3.17) 



38 



otj = common-base current gain in the inverse mode 

I~ 0 = saturation current of collector junction with zero 

emitter current 

I EQ = saturation current of emitter junction with zero 
col lector current 

b. Transistor Model 

Figs. 3.2a and 3.3a are representations of (3.16) and (3.17) 
for a n-p-n and p-n-p transistor respectively. Figs. 3.2b and 3.3b repre- 
sent the corresponding equivalent circuits if each diode is modeled as 
outlined in the preceding section. 

The reverse saturation currents are as follows: 



I 



CS 



I 



ES 



a I I C0 
1 " a I a N 

Veo 



(3.18a) 

(3.18b) 



If the positive senses of the currents are as indicated in 



Figs. 3.2 and 3.3, I^q and I^q are positive quantities for n-p-n transistors 
and negative quantities for p-n-p transistors and C<^(v) = 



and C $E (v) Tp E 



di 



dt 



'DC dt 



de . Tq E is approximated by the transistor rise 



time and T^ is approximated by the sum of the transistor storage time and 
the transistor fall time. The rise, fall and storage times are complex 
functions of the base current, collector current, a^, , and the transistor 

current-gain cut-off frequency, and are defined in Fig. 3.4 [17]. 

If C^^(v) = C j-jQ + C^q(v) and Cj^(v) - C n p + p ( v ) the 

equations for Fig. 3.2b become 



'DE °SE' 



Cj„(v) 



dv~ 



dc' ' dt 



+ 1 , + 7T 

dc R 



SHC 



Vde 



(3.19a) 
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(a) 



N de 



a i 'dc 




(b) 



V C 



* 



+ V E 



Fig. 3.2. Model and equivalent circuit for the n-p-n 
transistor. 
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®N 1 BE “i'bC 




(a ) 



a 



N 



de 



Vdc 




Fig. 3.3. Model and equivalent circuit for the p-n-p 
transistor. 
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Time 



(b) 



Fig. 3.4. Transistor switching times. 
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0 



(3.19b) 



dv E 

C de (v) dF 



'de 



SHE 



Vdc 



which are in the form of (3.3) where 

i 






f(x(t)) = 



a N 1 de 



a lV 



'dc 



'de 



R SHC C dc (v) 



( V ) ^CUF^rla( 



SHE de 



wt) 



(3.20a) 



B = [0] (3.20b) 

Similar equations can be written for Fig. (3.3b) 

The transistor example (III. E.) does not include the effects 
of parameter temperature dependence and C^, C^, and are assumed 

to be constant. Typical parameter values are used. 



D. DIODE EXAMPLE 



A diode problem of considerable interest is now solved in order to 
illustrate the modeling and the mathematical formulation discussed in this 
chapter. The problem considers the effect of a radiation current pulse 
injected across a p-n junction as indicated in Fig. 3.5, where the radiation 
current pulse is defined by 

di (t) 

ipp^) + t r dt = A (3.21) 



where y(t) is the unit pulse given in Fig. 3.6. 

Replacing the diode in Fig. 3.5 by the model described above yields 
Fig. 3.7. Typical values are 

1^ = 1.0 pA C D = 1.0 pF T Q = l.Ops 

a = 1/ .026 V -1 R sh = 100.0 M T R = 0.2ps 

A = 10.0 mA 
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Fig. 3. 5. Diode circuit. 



7(t) 



' time {ft s) 

Fig. 3.6. Radiation current 
pulse. 



R 



R 



SH 



i 



5 ^ 


C D 


Cc(v) ^ 

- c 


'd 







■pp 



Fig. 3.7. Equivalent circuit. 



44 



1 . Analytical Solution for Radiation Current 

Equation (3.21) can be solved using Laplace transforms. 



A* ( 



1 



1 



I (s) 

PP 



S + 1/Tr 



) 



0 < t < 10 
-7 



-7 



(3.22) 



• 394 A '< s - i/t r ) t i 10 



The resulting time response is given in (3.23). 



v (t) 



f A- (1 - e" t/T R) 



,394 A e -<t-10- 7 )/T R 



0 < t < 10 



t > 10 



-7 



(3.23) 



2. Piecewise-Linear Diode Model 

If the diode current- voltage relationship is assumed to be piecewise 

linear (Fig. 3.8) i d = 0 for v< 0 and v = i d R d for v > 0, where R d = R d (v) 

= r d (v). The storage capacitance becomes a constant, C $ = ^D/R d , for 
v > 0 and an open circuit for v < 0. Fig. 3.7 can be redrawn as in Fig. 

3.9 where Ry is the parallel combination of R, R $H and R d and C = C $ for 
v > 0, and R^ is the parallel combination of R and R^ and C = C D for v < 0. 

The rise and fall time can be determined using Laplace transforms. 



For 0 <t < 10 



Z(s) = 



1/C 



s + 1/RC 



Ks) = I pp (s) - 



E/R 

s 



(3.24a) 



(3.24b) 



= A.( 



1 



s + 1 /T e 



E/R 



The voltage takes the form 



v(t) 



C ] + C 2 e' t/T R + C 3 e" t/RC 



0 < t < 10" 



(3.25) 
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Fig. 3.8. Piecewise linear model . 




Fig. 3.9. Diode circuit representation . 
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For t > 10 



-7 



■s ■ • 394 A -<ft 1 vt 



R 



E/R 

s 



(3.26) 



The resulting voltage response is given by 



v ( t) = C 4 + C 5 e" t/T R + C 6 e" t/RC t >_ 10 



-7 



(3.27) 



The differential equation becomes 



= i - Gv - E/R 
dt pp 



(3.28) 



where G = 1/Ry. 

Equation (3.28) is solved using trapezoidal integration. Several 
values of R^ are used and the output voltage, v q = -v, is indicated in 
Fig. 3.10. The computer program is presented in Appendix B. The solution 
in Fig. 3.10 is of little practical value because of the wide variation of 
the voltage response with assumed values of R^. 

3. Taylor Series Approximation 

Equation (3.28) may be rewritten as 



C ( v ) dv = ipp dt - 



G( v) • v dt - E/R dt 



(3.29) 



where C(v) = Cp + C<.(v) 



Gtv) ' R + 4 + 



Using trapezoidal integration for (3.29) yields 



[C(v(n)) + C(v(n-1 ))] ^ ^ [1 (n) + i (n-1) - ^ ] 



at 



2E 



PP 



PP 



(3.30) 

- ^ [G(v(n) ) • v(n) + G(v(n-1 ) ) *v(n-l )] 
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OUTPUT VOLTAGE (volts) 




TIME (/is) 



FIG. 3.10. NONLINEAR EXAMPLE — PIECEWISE LINEAR 
MODEL WITH R d ADJUSTED FROM 
0.1 kXi TO 4.0 kft. 
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Expanding C(v(n)) and G(v(n))-v(n) in a Taylor series and including 
only the first-order terms yields 

C( v(n) ) = C(v(n-1)) + |£| • Av = C(v(h-1 ) ) (3.31a) 

av v(n-l) 

G(v( n) ) • v(n) = G(v(n-1 ))• v(n-l ) + G(v(n-l))*Av + |^-| ^Av«v(n-1 ) 

= G( v ( n-1 ) )• v(n) 

where it is assumed =0 and | = 0 . 

8v v(n-l) 8v v ( n-1 ) 

Substituting (3.31) into (3.30) results in 
C(v(n-l)).Av . ^[i pp (n) ♦ yn-D -f] 

- ^ [G( v{ n-1 ) ) • (v(n) - v ( n-1 ) ) ] 

Solving for v(n) yields 

v (n) = 2C(v(n-1)) AT v G(v(n-l)) v(n . 1} 

2C(v(n-l )) + AT*G(v(n-l )) 

i ( n ) + i (n-1 ) - -^f- 
+ aT . . P E B 

2C(v(n-l ) ) + AT * G(v(n-1 )) 

It is of interest to note that the same result is obtained if 
(3.28) is integrated over the interval (n-l)T to nT as a linear, constant- 
coefficient equation with C(v) and G(v) having the values C(v(n-1)) and 
G(v(n-1)) respectively. 

This method is implemented with an adjustable integration step 
size. The step size is adjusted to insure that C(v) and G(v) do not vary 
by more than a fixed percentage from (n-l)T to nT. The resulting output 
voltage, v o , is indicated in Fig. 3.11, where the maximum that C(v) and 



(3.31b) 



(3.32) 



(3.33) 
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10 . -• 




TIME (^s) 

FIG. 3 . 11 . NONLINEAR EXAMPLE -- DIODE RESISTANCE AND 
CAPACITANCE ARE HELD CONSTANT BETWEEN 

CALCULATION POINTS. PERCENTAGE VARIATION 
INDICATES THE MAXIMUM ALLOWED VARIATION 

IN C (v) AND G(v) BETWEEN CALCULATION POINTS. 
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G(v) are allowed to vary is 10%, 1% and 0.1%. The computer program is 
presented in Appendix B. This method becomes more accurate as the percent- 
age variation allowed for C(v) and G(v) is reduced. 



4. Diode Solution using General Nonlinear State-Variable Format 

The circuit in Fig. 3.7 contains a single state variable, v(t). 

The state equation for the circuit can be represented in the form of (3.3) 

where w /.n 

- G-v(t) - I e ( e av ^ _ i) 

(3.34a) 



f(v(t)) 

B(v(t) ) 

u(t) 
1 



C D + T D I S 0te 
1 



S__ 

av(t) 



C D + T D I S ae 

i (t) - E/R 
PP 



av(t) 



(3.34b) 

(3.34c) 



where G = ^ 



SH 

For the single state-variable problem considered in this example 
(3.9) reduces to 



v(u) 



AT [f + B (u(n) + u(n-l))] 

AT r 9f 3(B ♦ u(n)) -■ 

2 l 3v 3v' J 



(3.35) 



where the partial derivatives may be obtained from (3.34) and f, B, ■*— , 

dV 

and a v U ’ ^ n ^ are evaluated at v(n-l). 

The computer program for this method appears in Appendix B and the 
results are indicated in Fig. 3.12. The program adjusts AT automatically 
so that Av never exceeds 0.01 volts. These results have been confirmed 
by direct integration of the nonlinear state equation using a fourth-order 
Runge-Kutta integration scheme. 
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OUTPUT VOLTAGE (volts) 




TIME ( jj.s ) 



FIG. 3.12. NONLINEAR EXAMPLE — SOLUTION USING 
TAYLOR SERIES METHOD. 
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5. Presentation of Data 



The results of several runs are presented in Figs. 3.13 through 
3.17. The recovery time, defined as the time for the voltage to recover 
to 90% of its original value, is plotted as the ordinate. The abscissa 
represents the parameter value. 

In Fig. 3.13 the amplitude of the radiation pulse. A, is varied 
from 0.265 mA to 20.0 mA. When A <_ 0.265mA the radiation current pulse 
does not cause the voltage to drop below the 90% value. 

In Figs. 3.14 through 3.17 Cp, R, T^, and T^ are adjusted with 
A = 5.0 mA and A = 10.0 mA respectively while the other circuit parameters 
are held constant at their nominal values. 

E. TRANSISTOR EXAMPLE 

The effect of a radiation current pulse injected across the base- 
emitter p-n junction of a transistor operating in the active region (Fig. 
3.18) is considered in this section. The radiation current pulse is 
defined by (3.21) where y(t) is the unit pulse given in Fig. 3.19. 

1 . Transistor Solution using General Nonlinear State-Variable Format 

Modeling the transistor in Fig. 3.18 by the equivalent circuit 



presented in Fig 


1 . 3.2 results in Fig. 3.20 where 




C SC^ 


av c 

t dc *cs ae 


(3.36a) 


C SE (v) 


aVr 

“ t de 1 es ae 


(3.36b) 


n de 


- ’es (e “ VE • ’> 


(3.36c) 


1 dc 


- 


(3. 36d) 


Typical 


parameter values used are 




a |\j = 1*0 


I<-Q = O.OlnA 


Tq£ = 1 .Onsec 
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( stV) auiij Xj®Aooay 
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•265 5. 10. 15. A ( m A ) 20. 

Fig. 3.13. Diode recovery time as a function of radiation current pulse amplitude. 



Recovery ti me (/x. s ) 



1.2 




I. 2. 3. 4 . 5. 

C 0 (pF) 

Fig. 3.14. Diode recovery time as a function of C D . 
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Recovery time ( fis) 




Fig. 3.15. Diode recovery time os a function of R. 
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Recovery ti me [/ ms ) 



1.2 




Fig. 3.16. Diode recovery time as a function of T R . 
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Recovery time ( /x s) 
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Fig. 3.17. Diode recovery time as a function of T D 
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Collector 




Fig. 3.18. Transistor circuit. 



r(t) 



1.5 

time (fis ) 



Fig. 3.19. Radiation current pulse. 
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Collector 




Fig. 3.20. Equivalent transistor circuit. 
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cij = 0.99 


I se = 0-1 pA 


Tpj, = 0.1 ys 


T^ = 0.4 nsec 


Cqe = i 


a = 38.46 V 


R $hc = 100 


C DC ’ 1 - 0pF 




r she = 100 MQ 


A = 0 . 1 mA 





The circuit equations are obtained by summing the currents at the 
collector and emitter nodes yielding 



C * V E V E 

C de^ cfiT + R^ + 1 de " a lV 1 pp " 1 R = 0 
dv v r \/ 

^dc^ dt + + n dc " °Vde + n R + R^" ~ ^ 



where 



= r (v cc + v c- V E> 
c 



Cd c ( v ) = + ^ 



SC' 



'DC 



(3.37a) 

(3.37b) 

(3.37c) 



C de (v) = C SE (v) +C DE 



The biasing resistors, R^ and R c> are calculated for a desired 

biasing voltage V , base=col 1 ector voltage v r , and base-emitter voltage 
c dv E dv„ 

v £ under steady-state conditions, by setting and i 

equal to zero in (3.37). The program used is presented in Appendix B, 

and the resulting values are R, = 0.35 and R = 4.5 kft for 1/ = 10.0 

3 b c cc 

V£ = 0.6 V and v^ = -3.4 V. 

Substituting (3.37c) into (3.37a) and (3.37b) and solving for 

• • • T 

— = ^- V E V C^ yields a result identical in form to (3.3) with 



f(x(t)) 



TE 



C de^ >R c 



V, 



C dc (v),R c 



TC 
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(3.38a) 



S l CS 






/ T / ' 

'“N E3 te 



'dc 



(e“ v C - 1) I Ft; (e av E - 1) 

' ~de<’'> 

I sc (e“ Vc - 1)' 

' “dc'W- 



aV[r 

Tv] 



i) 



B(x(t) ) = 












(3.38b) 



u(t) 



i (t) 
PP 



cc 



(3.38c) 



where G TF = i- + p— ^ — 
TE c K SHE 



and Gyp „ + _ + ^ 

10 K b K c R SHC 



The iterative solution equation is of the form of (3.9) where 
Vf T and V(B_u(n))^ are obtained by performing the indicated gradient opera- 
tion on (3.38) . 

The computer program for this problem appears in Appendix B, and 
the results are indicated in Fig. 3.21. The program adjusts AT automatically 
so that the norm, |e| , never exceeds 0.01 volts. 

|e| = /(^ V E ) 2 + ( AV C> 2 (3.39) 

These results have been confirmed by direct integration of the 
nonlinear state equations using a fourth-order Runge-Kutta integration 
scheme. 



62 




TIME (^s) 



FIG. 3.21. TRANSISTOR RESPONSE TO RADIATION 
CURRENT PULSE. 
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2. Presentation of Data 



The results of several runs are presented in Figs. 3.22 through 
3.25. The recovery time is plotted as the ordinate and the abscissa 
represents the parameter value. 

In Fig. 3.22 the amplitude of the radiation pulse, A, is varied 
from zero to 3.0mA. When A <- 2.0yA the radiation current pulse does not 
cause the voltage to drop below the 90% value. 

In Figs. 3.23 through 3.25 Tp C , T^ , C DC , and T R are varied 
respectively while the other circuit parameters are held constant at their 
nominal values. 
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Fig. 3.22. Transistor recovery time as a function of radiation current pulse amplitude. 




Fig. 3.23. Transistor recovery time as a- function of T DE and T DC 
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Fig. 3.24. Transistor recovery time as a function of C DE and C DC 
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IV. SENSITIVITY ANALYSIS OF NONLINEAR CIRCUITS 



A. INTRODUCTION 

Sensitivity analysis of nonlinear circuits and systems is a relatively 
new field. Recently S. R. Parker [4] has considered the sensitivity analysis 
of general nonlinear circuits utilizing auxiliary coupled networks. The 
sensitivity analysis of nonlinear circuits utilizing the concept of state 
space is covered in the next section, followed by a sensitivity analysis of 
the nonlinear diode circuit discussed in Chapter III. 



B. SENSITIVITY-AUGMENTED STATE EQUATION 

The state equation format adopted in Chapter III (repeated below) is 
convenient when considering circuit sensitivity. 

x(t) = A(x(t))-x(t) + c_( x_( t ) ) + B(_x(t) )• u(t) (4.1) 

x i ax. 

As in the linear case, the sensitivity function, S a = ■ — , 

a > a 111 a * 

J vJ 

may be obtained by performing the indicated partial differentiation on 
(4.1), yielding the sensitivity-augmented state equation. 
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isj. 
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CO 
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M sj 

— -J 




A c . 1 B c . 

J ' l 




£ 



(4.2) 



where j = 1,2, ...,r and x Sj , Ap, Bp, A gj , B Sj , u Sj are defined in Chapter 



II and 



*sj 



3C, 



ac. 



a In a, a In a. 
j j 



_£n_ 
3 In a. 

J 



(4.3) 



Separating the set of first-order partial differential equations for 
each sensitivity parameter in (4.2) yields 
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x c . = A x_. + c_. + A c . x + B u c . + B c . u 
— S J ~ S J ~ S J S J - -sj sj - 



for j = 1 ,2, 



Equation (4.4) may be rewritten as 



2< S j = + B ^M S j + B sj (x,x sj )u 

where j = 1,2, ... ,r and 



(4.4) 
.. ,r 

(4.5) 



f . (x,x . ) = Ax. + c. + A.x 

-sj - -sj -sj -sj sj - 



Applying trapezoidal integration to (4.5) yields 



AT 



4j( n ) ‘ = 2 ~ Cf s j(x(n),x si (n)) + f s1 (x(n-l ) ,x si (n-l ))] 



sj 



sj 



sj 



-sj 



+ [B(x(n))u S j(n) + B(x(n-1 ))u s j(n-l )] + j- [B SJ .(x(n) ,x s j(n) ) u(n) 



+ B SJ .(x(n-l),x s j(n-l)) u(n-l)] 



(4.6) 



The functions f s j(x_(n) ,x S j(n) ) and B $ j(x_(n) ,x $ j(n))u_(n) , where ^(n) 
is considered a constant, can be expanded in a Taylor series about the 
point (x(n-l ) ,x,. .(n-1 )) and approximated by the first two terms, as in 
Chapter III , yielding 

f Sj .(x(n) ,x Sj .(n)) = f Sj -(x(n-l),x S j(n-l)) 

+ (AX T v x )f S j(x(n-l ) ,x S j(n-l )) + (axJjV^ )f s j(x.(n-l ) ,x S j(n-l ) ) (4.7a) 

B SJ -(x.(n) ,x S j(n))u_(n) = B $ ..(x(n-1 ) ,x S j(n-l ))u.(n) 

+ (Ax T v x )[B S j(x(n-l),x s j(n-l))u_(n)] + (AxJjV x$ j)[B S j(x(n-l ) , x s j ( n-1 ) ) ujn)] 

(4.7b) 



where 



v = [— — 

X 1 3X-, 



3X 



3 -] T 



= [■ 



V xs j - 3*1 



3X 



SJ 



n . 
sj 
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B(_x(n))- u .(n) , where u .(n) is considered constant, can be expanded 

b J b J 

in a Taylor series about the point >^(n-l) and approximated by the first 
two terms yielding: 



B ( xj n ) ) • u, ( n ) = B(x(n-l))*u .(n) + (Ax T V )[ B(x(n-1 ) ). u .(n) ] 

~ J — bj * bj 



(4.7c) 



Substituting (4.7) into (4.6) yields 



4j<"> + + 



-SJ 



>J XSJ'-SJ' 



+ [B^Ugjfn) + Usj(n-l)) + (Ax T V x ) (Bu^n) )] 



+ T [ B s j’(M.( n ) + H.(n-D) + (M Tv x )( B s jU.( n )) + (^^^(B^r,))] 

(4.8) 

where f B . and B .u(n) are evaluated at the point (x(n-l)), 

— sj sj sj— r — ' ’ 

Recall that 






(4.9) 



Substituting (4.9) into (4.8) and solving for x.(n) yields 

^ J 

%<"> - x sj ("-l)*AT[I-^ t (V XSJ f s T J ) T + (V xs /B 5 .u(n)f) T ) 



[f 4 + i {(Ax T V )f . + B* (u (n) + U (n-1 ) ) + (Ax T V ) [Bu (n)] 

“SJ C- x “SJ “SJ “SJ X “SJ 

+ B *(u(n) + U_( n - 1 ) ) + (Ax T V )[B _u_(n)]} ] (4.10) 

S J A b J 

Equations (3.9) and (4.10) can be solved directly for _x(n) and x.(n) 

5 J 

from the known values of x_(n-l), ujn), ujn-1), u .(n), and u .(n-1). The 

b J b J 

solution yields the time response of the states and their corresponding 
sensitivities. 
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C. NONLINEAR EXAMPLE 



In order to illustrate the technique discussed above, sensitivity 
analysis of the diode circuit in Chapter III is considered. The state 
equation for the diode circuit of Fig. 3.5 contains a single state var- 
iable. Equation (4.10) reduces to 

v s j (n ) = V sj(n-D + iT -[ f so + iv + B-(u sj (t>) + u sj (n-D) 



3(B*u .(n) ) 



sj 



3(B .u(n)) 



3v 



Av + B s j‘( u ^ n ^ + L, ( n - 1 )) + — 



Av}] 



n AT , 3f sj . 3(B sj u < n >> , n -l 

u ~ 2 *3v 



sj 



3v 



} ]' 



sj 



where f 



df. 3(B*u (n)) 3(B _.u(n)) 

sj* V • B ’ B sj> §v and — V— 



(4.11) 



are evaluated 



at the point v(n-l ) . 



The state equation may be written in the form of (4.1) where 



A(v(t) ) 
c (v(t)) 



C D + I D I s ae aV 
I s (e av - 1) 

C D + Vs 0 *™ 



(4.12a) 

(4.12b) 



and B( v(t) ) and u(t) are given in (3.34b) and (3.34c) respectively. 

The equations for the sensitivity functions take the form of (4.5). 
The sensitivity parameters considered are the diode parameters, Tp and 
Cp, the radiation-pulse parameters, A and T^, and the circuit parameters, 
R and E. The terms in (4.5) are presented in table 4.1 for each of the 
sensitivity parameters where 

k l = I S ae ° t k 2 = T D k l 
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Table 4.1. Terms of (4.5) for the circuit of Fig. 3.6. 



^3 = a ^2 



k 4 = I s (e av -1) 



D 



C D + k 2 



3i 3i 

The expressions for — ^ — and y y -P P - y are obtained by performing 

3, " A R 3f. ,(B.u s1 (n» 

the indicated partial differentiation on (3.21) and — - J , ^ 



3fB s,i u(n)) 3f s,i 



3v 



3V Sj 



, and 



8(B sj u(n)) 

3v~ 

sj 



are obtained by performing the 



indicated partial differentiation on the equations listed in table 4.1- 
In addition to the sensitivity functions the voltage response, 
v = -v , is shown for (1) nominal parameter values, and (2) a 10% change 
in the sensitivity parameter value. The response is presented in Figs. 
4.2, 4.4, 4.6, 4.8, 4.10, and 4.12. The computer program for generating 
the sensitivity functions is given in Appendix B, and the resulting 
sensitivity functions are presented in Figs. 4.1, 4.3, 4.5, 4.7, 4.9, and 
4.11. The program adjusts AT automatically so that |c.| never exceeds 

J 

0.05 volts and I Av | never exceeds 0.01 volts, where 



(Av) 2 + (Av sj ) 2 



for j = 1 , 2, ... ,6 



(4.13) 



These results have been confirmed by direct integration of the non- 
linear state equations using a fourth-order Runge-Kutta integration scheme. 
It should be noted that the method outlined in this chapter is approxi- 
mately twice as fast as the fourth-order Runge-Kutta integration scheme. 



D. INTERPRETATION AND APPLICATION OF THE SENSITIVITY FUNCTION 

The definition of the sensitivity function adopted in this thesis is 
x. 3x . 3x. 

S a ; = 9 i n n a = ^ a ./ a j . As noted, this definition applies only for 

J J J J 

infinitesimally small changes in x^ and a^. The sensitivity function may 
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FIG. 4 . 1 . SENSITIVITY FUNCTION WITH RESPECT TO T D . 
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FIG. 4.2. RESPONSE OF DIODE CIRCUIT TO RADIATION 
PULSE. 
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FIG. 4.3. SENSITIVITY FUNCTION WITH RESPECT TO E. 
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FIG. 4.4. RESPONSE OF DIODE CIRCUIT TO RADIATION 
PULSE. 
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FIG. 4.5. SENSITIVITY FUNCTION WITH RESPECT TO R* 
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FIG. 4.6. RESPONSE OF DIODE CIRCUIT TO RADIATION 
PULSE. 
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FIG. 4.7. SENSITIVITY FUNCTION WITH RESPECT 
TO A. 
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FIG. 4.8. RESPONSE OF DIODE CIRCUIT TO RADIATION 
PULSE. 
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FIG. 4.9. SENSITIVITY FUNCTION WITH RESPECT TO Cd. 
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FIG. 4.10. RESPONSE OF DIODE CIRCUIT TO RADIATION 
PULSE. 
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FIG. 4 . 11 . SENSITIVITY FUNCTION WITH RESPECT TO T R . 
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FIG. 4 . 12 . RESPONSE OF DIODE CIRCUIT TO RADIATION 
PULSE. 
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be used to approximate changes in a response due to fractional changes 
in a parameter as follows 



AX i 




•( 




(4.14) 



In general this approximation becomes more accurate as the percentage 
change in the parameter value approaches zero. Table 4.2 represents a 
comparison of results obtained for the underdamped response of the linear 
example in Chapter II. The output voltages, v 0 » v ] * ••• , v 5 , are 
obtained using resistance values of R Q , 1.1 R q , ... , 1.5 R q , where R q 
is the nominal value. Approximations for these voltages are calculated 
using (4.14); specifically, 



w AR „ 

v a v o + S R '< TT 1 < 4 - 15 > 

0 

where v fl is the approximate response after R q has been changed by AR Q , 
and v Q is the nominal response. Table 4.3 presents a similar comparison 
for the nonlinear diode example. Values are tabulated for time intervals 
when the sensitivity function does not vary rapidly. The results in 
tables 4.2 and 4.3 indicate that the approximation for the change in output 
voltage improves as the amplitude of the sensitivity function approaches a 
constant value. In this region correct values are predicted for variations 
as large as 50%. It is noted that (4.14) does not yield useful results 
for parameter variations as small as 1% when the sensitivity function dis- 
plays sharp peaks. In this region the sensitivity function may still be 
used to determine changes in circuit responses by assuming that the asso- 
ciated parameter variations cause a time shift rather than a voltage 
change. Fig. 4.13a represents an idealized switching-circuit response 
where the effect of an incremental parameter variation is upon the switch- 
ing time. Fig. 4.13b represents the corresponding sensitivity function. 
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Table 4.2. Voltage comparison for linear example using 
(1) actual parameter changes, (2) approximation using (4.15). 
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Table 4.3. Voltage comparison for nonlinear diode 
example using (1) actual parameter changes, (2) approximation 

using (4.15). 
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FIG. 4.13 (a) IDEALIZED CIRCUIT RESPONSE 
AND (b) ASSOCIATED SENSITIVITY 
FUNCTION. 
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The delay time, T^, due to the parameter variation is given as 




(4.16 y 



Applying (4.14) to Fig. 4.13 it follows that 




(4.17) 



The area, A, under the sensitivity function over the time interval 
t 2 - t-j is 

t x 

A = / S. 1 dt = h T . (4. 

t, a j a 

Solving for h in (4.17) and substituting into (4.18) yields an 
equation for the delay time T^. 



Equation (4.19) may be applied to estimate the delay time correspond- 
ing to a pulse-type sensitivity function by calculating the area of the 
pulse and substituting into (4.19). This type of calculation is performed 
for the waveforms of Figures 4.1, 4.3, 4.5, and 4.7 for 10%, 20% and 50% 

parameter variations. The results are compared with the actual time 

t 2 x. 

delays in table 4,4. The area, A = / S dt , is calculated assuming 

ti a j 

that the amplitude spikes are triangular. Times t, and t ? are chosen at 

x. 

those points where S 1 = .15 h . The results indicate time delays can be 

j 

calculated within 10% for parameter variation as large as 50%. It is 
noted that calculation of the area presents the greatest limitation in the 
accuracy of the delay time. 



T 



d 




(4.19) 
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times for given percentage changes in parameter values. 



V. CONCLUSIONS 



The technique of generating sensitivity functions for linear circuits 
using the concept of a sensitivity model is shown to be an extension of 
the compensation theorem. A set of iterative equations is developed to 
yield simultaneously the time responses of the circuit states and of their 
corresponding sensitivities. An example illustrates the procedure. 

Diode and transistor modeling for discrete computation based upon 
the Ebers-Moll equations is considered. The models adopted contain both 
nonlinear capacitors and nonlinear resistors for the junctions. A state- 
variable formulation is developed which results in a set of nonlinear 
iterative equations for circuits containing linear elements, diodes, and 
transistors. The use of these equations is illustrated by considering the 
effects of a radiation pulse on a p-n junction, and the resulting recovery 
times of the junction are tabulated. 

The foregoing nonlinear iterative equations are extended to generate 
simultaneously the sensitivity functions as well as the circuit solution. 

The preceding example, which considers the effects of a radiation pulse on 
a p-n junction, is then analyzed from the point of view of a sensitivity 
analysis. This analysis indicates that proper interpretation of the sen- 
sitivity function results in accurate estimated for changes in circuit 
responses due to parameter variations as large as 50%. When the amplitude 
of the sensitivity function is constant the change in the circuit response 
due to a parameter variation is reflected as a change in response ampli- 
tude. However, in regions where the sensitivity function displays impulse- 
like amplitude peaks, the change in the circuit response due to a parameter 
variation is reflected as a time delay or shift in the response characteristic. 
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APPENDIX A 



PROOF OF (ax T V) f = (Vf T ) T AX 



The matrix equations Ax Vf may be written as 



(Ax T V) f = [ Ax-| Trrr- + Ax„ rX- + 



3x, 



2 3x, 



AX 



AX 



3fn 



1 3X-| 
3f, 



1 3x n 



+ AX 



n 3x„ 





3f > ♦ 

3x 2 




3f 


AX 2 


' • + Ax n 


3X, 




3fp 




3f. 


Ax 2 


3x 2 + • 


. . f Ax n 





n 

2 I 



n 



1 



AX 



3f_ 



1 3x n 



3f 



+ AX 



2 3x , 



3f 



+ AX 



n 3x„ 



(A. 1 ) 



Consider this matrix equation 



Vf 



3X n 



[fl 



3X„ 






3 f l 


3f 2 


3f 

... n 


3X i 


3X-| 


3x-| 


3 f l 


3f 2 


3f 

. . . n 


3*2 


3X 2 

C 


3X 2 


3 f l 


3f 2 


... 3f n 


^ X n 


° n 


3x n 



(A. 2) 
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Taking the transpose of equation (A. 2) and multiplying by Ax yields 



(Vf T ) T Ax 



k 


3f l 




» f i’ 


i 

1 


j 


3X-| 


3X2 




3x n ; 




j Ax ] 


3f2 


3f2 




3f 2 




! Ax 2 




3X2 




n 




• 


• 




• 




. 


• 


• 




• 




* 


3f n 

3x l 


9f n 

3X2 




3f n 

8x n 






Q • 0 




Ax 

n_ 



r 



Ax 



9fj 
1 3X n 



AX 



3fj 

2 3X. 



AX 



ih 

n 3x„ 



Ax 



1 3X n 



AX 



% 

2 3X. 



Ax 



3 f 2 
n 3x~ 



Ax 



!!r 

1 3X, 



AX 



Hr 

2 3X, 



AX 



5 

n 3x„ 



Equating equation(A.3) to equation(A.l) results in equation (A. 4), 
completes the proof. 

(Ax T V) f s (Vf T ) T Ax 



(A. 3) 



wh i ch 



(A. 4) 
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APPENDIX B 



COMPUTER PROGRAMS 



1. SOLUTION AND SENSITIVITY ANALYSIS PROGRAM FOR LINEAR 
CIRCUITS 



C 

c 



c 

c 



c 

c 

c 

c 



R F AL *8 ITITLE 

DIMENSION AI(2,2),AS(2,2),AMI(2,2),API(2,2), 

1AM I INV<2,2 ) ,PHI ( 2 ,2) , AMIT2( 2 ,2 ) ,BS ( 2, 1 ) , 

2 G AM ( 2 » l),BSl(2,l),BS2(2,l),BS3(2,l), 

3GAM2 1(2,1) »G AM22 ( 2 , 1 ) , GAM2 3 < 2 , 1 ) , G AM 1 1 ( 2 , 2 ) , 

AG AM 1 2(2,1) »GAM 13(2,2) , TEMP 1 ( 2 , 2 ) , TEMP 2 ( 2 , 2 ) , 
5TFMP3(2, 2) ,TEMP4( 2,2) ,US< 1, 1 ) ,XS< 2, 1 ),XS1< 2, 1) , 
6XS2( 2, 1) ,XS3( 2 ,1 ) , TEMPS1 (2, 1 ) ,TEMPS2(2 ,1 ) , 
7TFMPS3(2, 1 ) ,TEMPXS( 2, 1 ) , S 1 < 910) ,S2 (910 ) , S3 (910) , 
8XP(91C),YP(910),AS1(2,2),AS2(2,2),AS3(2,2), 
9ITITLE (12) 

T I ME = C .C 
US MATRIX = 2*U 
US ( 1, 1 )= 14C. 

EQUATIONS OF THE FORM 
X ( NT ) = PHI*X((N-1)T) 

XSJ(NT) = PHI*XSJ((N-1)T) 

+ GAM 1 J* ( X ( NT ) +X ( (M-l)T) 

WHERF 

AMI , 



+ GAM* ( U ( NT ) +U ( (N- 1 ) T ) ) 

+ GAM*(USJ(NT)+USJ( ( N — 1 ) T ) ) 
+ GAM2J*< U< NT )+U( (N-l) T) ) 



( I-(T/2)*A) 



-1 



(I+(T/2)*A) = API 



(I-(T/2)*A) = AMIINV 

-1 



PHI 



( I + ( T /2 ) A ) 



<I-(T/2)*A) (I+(T/2)*A) 

-1 

( I- ( T/2 ) *A ) (B)*(T/2) = GAM 

-1 

( I - ( T /2 ) * A ) ( A S J* ( T/2 ) = G AM1J 

-1 

( I — ( T/2 ) * A) (BSJ)*(T/2) = GAM2J 

NS = U OF STATES 
NU * # OF FORCING FUNCTIONS 

FORMAT AND DIMENSION STATEMENTS MUST BE CHANGED WHEN 
THE NUMBER OF STATES OR FORCING FUNCTIONS ARE CHANGED 
NS=2 
NU= 1 

R FAD IN THE I AND A MATRICES 

READ (5,1) ((AI(I,J),J=1,NS),I=1,NS) 

R E AD ( 5 , 1 ) ( ( A S ( I , J ) ,J = 1,NS) ,1=1, NS) 

1 FORMAT (8E1°.2) 

CALCULATE 

AM I = ( I-( T/2) A) API 

T=. 16E-03 
Q = T/ 2 . 

CALL CONST (Q, AS, NS, NS, TEMPI) 

CALL SUB(A I, TEMPI, NS, NS, AMI ) 

CALL ADD(A I, TEMPI , NS, NS, API ) 

CALCULATE -1 

AMIINV = ( I — (T/2 ) A ) 

CALL REC IP ( .000001 , AM I , AM I I NV , KE R , NS ) 

I F (K EF -2 ) 2,3,2 

3 WRITF(6,4) 

4 FORMAT (5HKER=2 ) 

GO TO 4C 

2 CONTINUE 

CALCULATE PHI = AMIINV*API 
CALL PROD (AMIINV, API, NS, NS, NS, PH I ) 

CALCULATE AMI T2 = AMIINV*(T/2) 

CALL CONST (Q, A MI I NV, NS, NS , AMI T2 ) 

READ IN B MATRIX 

RE AD (5,1) ((BS(IfJ) » J = 1 » N U ) ,1=1, NS) 

CALCULATE GAM = AMI T 2 *B 
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CALL PROD ( AMIT2»BS»NS*NS»NU»GAM) 

C READ IN BS MATRICES AND CALCULATE GAM2J = AMIT2*BSJ 

RE AD (5,1) ( ( B S 1 ( I*J)»J=l»NU)»I=l,NS) 

RE AD (5,1) ( (BS2( I »J) , J = 1 , NU ) , 1 = 1, NS) 

Rf AD( 5 , 1 ) ( ( B S 3 ( I , J ) » J = 1 ,NU ) ,1 = 1, NS) 

CALL P ROD( AM IT 2, B S 1 , NS , NS , NU , GAM2 1 ) 

CALL P ROD ( AMIT2»BS2»NS»NS»NU,GAM22 ) 

CALL PROD ( AMIT 2, BS3»NS»NS»NU,GAM23) 

C READ IN AS J MATRICES AND CALCULATE GAM 1 J = AMIT2*ASJ 

READ (5,1) ((AS1( I ,J) ,J=1,NS) ,1=1, NS) 

RE AD (5,1) ((AS2(I,J),J=l,NS)tI=l,NS) 

RE AD (5,1) ( ( A S3 ( I ,J) , J = 1 , NS ) ,1 = 1, NS) 

CALL PROD( AMIT2,ASI,NS,NS,NS,GAM11) 

CALL P ROD ( AMIT 2, AS2 , NS , NS , NS , GAM 1 2 ) 

CALL P ROD ( AMIT2, AS3 ,NS ,NS , NS, GAM1 3 > 

WRITE! 6, 12 ) 

12 FORMAT ( 1 1* ,21X ,• A MATRIX') 

WRITF( 6, 14) ( ( AS( I , J) , J = 1 ,NS) , 1=1 , NS ) 

14 FORMAT (//, 2( 10X.E10. 2),/) 

WRITE(6,15 ) 

15 FORMAT ( //2 2X , ' B MATRIX') 

WRITE(6, 16 ) ( ( BS ( I * J ) * J= 1 ,NU I , 1 = 1, NS) 

16 FORMAT (//,1(21X,E10.2,/) ) 

WP ITF ( 6, 13 ) 

13 FORMAT (//22X, • PHI MATRIX') 

WPITF ( 6, 14 ) ( ( PHI ( I , J ) , J=1 ,NS) , 1 = 1 ,NS) 

WRITE(6»24) 

24 FORMAT (//22X, ' GAM MATRIX') 

WPITE( 6,16) ( ( GAM ( I , J ) , J = 1 ,NU ) , I = 1 , NS ) 

WR IT E ( 6 , 17 ) 

17 FORMAT!//, T27, 'SENSITIVITY MATRICES', /,T19,'GA MU' , T49 
1, • GAM2 1 ' ) 

WR ITF (6, 19) (( ( GAM 1 1 ( I »J) »J=1,NS) , (GAM21( I , J ) , J= 1 , NU ) ) 
1, 1=1, NS) 

19 FORMAT (//, 2(5X,E10.2,5X),5X,E10.2,/) 

WRITF(6,20 ) 

20 FORMAT!/// , T 19 . • G A M12 • , T49 , • G AM2 2 • ) 

WR ITE ( 6, 19 ) (( (GAM 12 ( I , J ) , J= 1 , N S ) , ( GAM 2 2 ( I , J ) , J= 1 , NU ) ) 
1,1=1 , NS) 

WRITE( 6,22 ) 

2 2 FORMAT (/// ,T 19, 'GAM13' ,T49, ' GAM2 3 ' ) 

WRITE (6, 19 ) ( ( ( GAM13 ( I , J ) , J= 1 , NS ) , ( GAM23 ( I , J ) , J = 1 , NU ) ) 
1, 1=1, NS) 

WR ITE ( 6, 31 ) 

31 FORMAT ( ' 1' ,/ / ,T13, 'OUTPUT' ,T 30 , • SE NS I T I V I T Y • , T 50 , ' SENS 
1 IT IVITY' ,T70,' SENSITIVITY' ,/,T33, 'WRT R',T53,'WRT L',T 
273, 'WRT C' ,T93,'TIME' ,///) 

C INITIALIZE 

RE AD (5, 26) ( ( X S ( I , 1 ) , I = 1 , NS ) , ( XS 1 ( I , 1 ) , I = 1 ,N S ) , ( X S2 ( I , 

11) ,1=1 ,NS) ,(XS3(I ,1) ,1=1, NS) ) 

RE AD (5,26) XP ( 1 ) , YP ( 1 ) , S 1 ( 1 ) , S2 ( 1 ) , S3 ( 1 ) 

26 FORMAT (3E1C. 2 ) 

1 = 1 
Q= 1 . 

30 WPITE(6,32 ) XS ( 1 , 1 ) , XS 1 ( 1 , 1 ) , XS2 ( 1 , 1 ) , XS3 ( 1, 1 ) , T I ME 

32 FORMAT! 5F20. 5) 

CALL CONST ( Q , XS, NS, NU,TEMPXS) 

CALL CONST (Q,XS1»NS»NU,TEMPS1) 

CALL CONST(Q,XS2,NS,NU,TEMPS2) 

CALL CONST (Q,XS3»NS»NU,TEMPS3) 

C CALCULATE STATE VECTOR 

CALL PROD( PHI , TEMP XS, NS, NS, NU, TEMPI) 

CALL PROD( GAM, US, NS, NU , NU , TEMP2 ) 

CALL ADD(TEMP1,TEMP2,NS,NU,XS) 

CALL ADD(XS,TEMPXS,NS,NU,TEMP4> 

C CALCULATE SENSITIVITY VECTORS 

CALL PROD ( PHI , TEMP SI , N S , NS , NU , TE MP 1 ) 

CALL PROD( GAM1 1 , T EMP4 , NS , N S , NU, TEMP2 ) 

CALL PROD ( GAM2 1 , US , NS , NU , NU , T E MP3 ) 

CALL ADD (TEMPI ,TEMP2, NS, NU, TEMPI) 

CALL ADD(T EMP1 »TEMP3»NS, NU , X S 1 ) 

CALL PROD ( PHI , TEMPS2 , NS, NS,NU, TEMPI ) 
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c 



CALL PP0D( GAM 12, T E MP4 , NS , N S , NU, TEMP2 ) 
CALL P ROD ( GAM2 2,US ,NS, NU, NU , TEMP3 ) 
CALL ADD (TEMPI ,TEMP2, NS, NU, TEMPI ) 

CALL A DO (TEMPI ,TEMP3,NS,NU»XS2) 

CALL PRCID( PHI , TEMP S3, NS, NS, NU, TEMPI) 
CALL PROD( GAM13,TEMP4,NS,NS,NU,TEMP2) 
CALL P ROD ( GAM2 3 , US , NS , NU , NU, TE MP 3 ) 
CALL ADD (TEMPI ,TEMP2, NS, NU, TEMPI ) 

CALL ADD (TEMPI »TEMP3»NS»NU,XS3) 



1 = 1 + 1 

TI ME = T I ME + T 
XP( I ) = T I ME 
YP ( I ) = XS ( 1 ,1 > 

S1(I) = XS1( 1,1) 

S 2 ( I > = XS2( 1, 1) 

S3 ( I ) = X S 3 ( 1,1) 

I F ( I .LT.9C 5) GO TO 30 
PLOT THE DATA POINTS 
RE AD ( 5 , 50 ) ITITLE 

50 FORMAT (6A8 ) 

RF AD ( 5 , 51 ) LAB ELI 

51 FORMAT ( A4 ) 

CALL DRAM 900, XP , YP , 0 , 0 , LABEL 1 , ITITLE, 3.E-C2, 10 

11.0. 2.2.5.6.1,L) 

WR IT E ( 6, 52 ) L 

52 FORMAT (//, IX, *L=' , 110) 

RE AD ( 5 , 5C ) ITITLE 

CALL DRAW! 900, XP, S3, 1 , C, L ABEL 1 , I T I TLE , 3 . E- 02 , 1 0 

12.0. 2.2.5.6. 1 , L ) 

WR I TE ( 6, 52 ) L 

CALL DRAW (900,XP,S2,2,0,LABEL1,ITITLE,3.E-02,10 

12.0. 2.2.5.6. 1 , L ) 

WP ITE ( 6, 52 ) L 

CALL DRAW ( 900, XP, SI, 3,0, LABEL 1, ITITLE, 3.E-C2,10 

12.0. 2.2.5.6.1.1) 

WP I TE ( fc, 52 ) L 

40 RETURN 
END 



• , 



• , 



• , 



• , 



SUBROUTINES 



SUBROUTINE ADD (A,B,N,M,C) 

C THIS SUBROUTINE ADDS TWO NXM MATRICES 

DIMENSION A(N,M) »3(N»M),C(N,M) 

DO 152 1 = 1, N 
DO 152 J=1,M 
152 C( I, J) =A( I ,J> +B( I , J) 

RETURN 

END 



SUBROUTINE SUB (A,B,N,M,C) 

C THIS SUBROUTINE SUBTRACTS TWO NXM MATRICES 

DIMENSION A(N»M) ,B(N,M),C(N,M) 

DO 152 1 = 1, N 
DC 152 J=1,M 
152 C( I , J)=A( I , J )- B( I , J) 

RETURN 

END 



SUBROUTINE CON ST ( Q, A , N, M, C ) 

C THIS SUBROUTINE MULTIPLIES A NXM MATRIX BY A CONSTANT 

DIMENSION A(N,M) ,C(N,M) 

15 DO 150 I = 1 , N 
00 150 J=1,M 
150 C( I, J)=Q*A( I , J ) 

RFTURN 

END 
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SUBROUTINE PROD ( A , 6 , N , M, L , C ) 

C THIS SUBROUTINE POST MULTIPLIES A NXM 

C BY A MXL MATRI X 

DIMENSION A(N.H),B(M,L)tC(N,L) 

DO 1 1=1, N 
DO 1 J=1 , L 
1 C( I, J)=C. 

DO 151 1=1 »N 
DO 151 J=1,L 
DO 151 K= 1 , M 

151 C( I ,J)=C( I ,J) + A(I , K ) * B ( K , J ) 

RETURN 

END 



C 



1 

2 

10 



1 1 
12 
13 



14 



15 

20 

30 

31 



32 

33 

35 

36 

34 
40 



41 

42 
4 3 

50 



SUBROUTINE REC I P ( E P , A , X , KER , N ) 

THIS SUBROUTINE TAKES THE INVERSE OF 
DIMENSION A(N,N) ,X(N,N) 

DO 1 1=1, N 
DO 1 J = 1 , N 
X ( I , J ) =C . 

DO 2 K= 1 , N 
X ( K,K) =1 . 

DO 34 L= 1 , N 
K P=0 
Z = 0. 

DO 12 K=L , N 

IF(Z.GE.ABS(A(K,L) )) GO TO 12 
Z=ABS(A(K,L) ) 

K P — K 

CONTINUE 

IFIL.GE.KP ) GO TO 20 
DO 14 J = L » N 
Z= A ( L , J ) 

A ( L, J ) =A ( KP, J ) 

A ( KP , J ) = Z 
DO 15 J= 1 , N 
Z = X ( L , J ) 

X ( L, J) =X(KP, J ) 

X(Kr J)=Z 

IF (ABS(A(L,L)) .LE.EP) GO TO 50 

IF (L.GE.N) GO TO 34 

L P 1=L + 1 

DO 36 K=L P 1 , N 

I F ( A ( K , L ) . EO.O. ) GO TO 36 

RATIO=A(K,L)/A(L,L) 

DC 33 J=L P 1 , N 

A(K, J)=A( K,J) -RATION AIL, J) 

DO 35 J= 1 , N 

XIK, J)=X(K,J)-RATIO*X (L, J) 

CONTINUE 
CONTINUE 
DC 43 1=1, N 
I I =N+ 1-1 
DO 43 J= 1 , N 
S = 0. 

IF ( I I .GE.N) GO TO 43 
I I P 1= I I + 1 
DO 42 K=I I PI ,N 
S= S+ A ( II ,K )*X( K. J ) 

X ( 1 1 , J )=( X I I I , J ) - S ) /A < I I , I I ) 

KER= 1 
RETURN 
K E R=2 
RETURN 
END 



MATRIX 



NXN MATRIX 
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DIODE SOLUTION USING PIECEWISE LINEAR MODEL 



C 



D I MENS ION V ( 60 0 ) ,CURR(600) , YP ( 905 ) , XP ( 90 5 ) 

REAL LABEL/4H / 

REAL*8 I T I TL E ( 121/' DIODE SOLUTION 

1 PIECEWISE LINEAR MODEL T=l.E-09 */ 

I C NT= 1 
I END=5 

RE AO (5, 3) SI ,E,CO,RSH,R,TR,TD,RADC 

3 FORMATI8EI0.2) 

100 T I ME=0 .0 

RF ADI 5 ,300 ) T , RD , NST OP, V ( 1 ) , XP (II , YP ( 1) , CURR ( 1 ) , TEMT 
300 FORMAT (2E10.2 ,1 10,5E10.2) 

ALP=L./.026 
GC LD=R ADC 
KKK= 590 
KKK I N= KKK- 1 
NM 1= 1 
6 N= 1 

4 IFITIME.LT. l.E-97) GO TO l 
GNEW=0 .0 

GO TO 2 

1 GNEW = R ADC 

2 N= N+ 1 



CURR ( N ) = ( ( 2 . *T R-T ) /I2.*TR+T) ) *C URR ( NM 1 ) + < T / ( 2. * TR + T ) ) * 
1 I GOLD+GNEW ) 

IFIVINMl) .LT.O.O) GO TO 9 
G=l./R+1 ./RSH+l./RD 
CAP=TD/P D 
GO TO 5 

9 G= 1 . /R +1 • / PSH 
C AP=CD 

5 V ( N ) = ( <2.*CAP-G*T)*V(NM1)+T*(CURR(N>+CURR( NM1 )-2.*E/R> 
1 )/ <2.*CAP+G*T ) 

GGLD=GNE W 
T I ME = T IME + T 

IFIN.GE.NSTOP) GO TO 20 
TEMT=TEMT+T 

IFITEMT.GF ..15E-C8) GO TO 20 
NM1=N 
GO TO A 
20 KKK=KKK+ 1 

I I =KKK-KKK IN 
XP ( I I ) =T I M E 
YP ( I I ) =-V ( N ) 

TE MT =0 . 

NM 1 = N 

IFIKKK.LT. 1A93 ) GO TO 6 
PLOT THE DATA POINTS 
I F ( IEMD.EO.l ) GO TO 104 
IF I ICNT.GT .1 ) GO TO 102 

CALL DRAW! 900, XP, YP,1 ,0, LABEL1 , IT ITLE , 3. E-07 , 3 . , 1,0,2, 
1C, 5, 5, 1 » L ) 

WRITEI6, 13) L 
13 FORMAT!//, IX, * L=» ,110) 

GO TO 1C3 

102 IF ( ICNT. EO. I END) GO TO 101 

CALL DRAW! 900 , XP , YP , 2 , 0 , LAB EL 1 , 1 T I TLE , 3 . E-07 , 3 . , 1 , 0 , 2 , 

10.5.5. 1 , L ) 

WRITE (6, 13) L 

1C3 ICNT=ICNT+1 
GO TO ICC 

101 CALL DRAW(900»XP»YP»3,0,LABEL1»ITITLE,3.E-07,3.,1,0»2, 
10,5,5, 1 , L ) 

WR IT E ( 6, 13 ) L 

GO TO 106 

104 CALL DRAW ( 900, XP»YP,0,0,LABEL1, ITITLE ,3. E-07, 3. , 1,0,2, 

10.5.5.1, LI 

WR I TE ( 6, 13 ) L 
106 RETURN 
END 
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3. DIODE SOLUTION HOLDING DIODE CAPACITANCE AND RESISTANCE 
CONSTANT BETWEEN CALCULATION POINTS 



RFAL*8 ITITLE 

DIMENSION V(600) ,CURR(60C) ,XP<905) ,YP(9C5) , 

I I T ITLE ( 12 ) 

I C NT= 1 
I E ND = 3 

RE AD ( 5 » 1 1 1 ITITLE 

11 FORMAT ( 6 A 8 ) 

RE AD ( 5 » 1 2 1 LAB EL 1 

12 FORMAT ( A4 1 

READ! 5,31 SI , E ,CD , RSH , R, TR , TO , RADC 
3 FORMAT (8F1C. 2 1 
100 T I MF=C .0 

CALL CANCEL! 2 ) 

READI5 ,3C0 ) T , GAM , NST OP 
300 FORMAT ( 2E 1 0. 2 ,110) 

TEMT=C .0 
ALP=1./.C26 

C INITIALIZE RD AND CAP 

P D= 1 . F +50 
C A P=CD 
V ( 1 1 =- 1C . 

CURR ( 1 )=C.O 
GOLD=R ADC 
KKK= 590 
KKKI N=KKK- 1 
NM 1= 1 
N= 1 

A IF(TIME.LT .l.E-07) GO TO 1 
GN EW=0 .0 
GO TO 2 

1 GNEW=R ADC 

2 N= N+ 1 
T=2.*T 

lF(T.GE.l.E-09)T=l.E-09 

9 CURR (N) = ( I 2 . *T R-T 1 / ( 2 . #TP+T 1 1 *CURR(NM1 1 + ( T / ( 2 . *TR +T ) 1* 
1 ( GOLD+GNEW 1 
G=1./R+1./PSH+1./RD 

V ( M = ( <2.*CAP-G*T)*V(NM1 1 +T* ( CURR ( N 1 +CUR R < NM1 1-2 .*E/R ) 
11/(2. *CAP+G*T 1 

RD=V(N)/(SI’MEXP(ALP*V(N) 1-1.1) 

GTTEST-l./R+l. /RSH+1. /RD 
CTTEST=CD+TD*S I*ALP*EXP( ALP*V(N) 1 
C0N1 = ABS ( G-GTT EST 1 
CCN2=ABS(CAP-CTTEST) 

CCN3=GAM*GTTEST 

C0N4=GAM*CTTEST 

IF( CONI. LT.CON3.AND.CON2.LT. CON A) GO TO 7 
IF(GAM.EQ.O.O) GO TO 7 
T= T/2 . 

GO TO 9 
7 GOLD=GNFW 
C A P=C T TE ST 
T I ME = T IME + T 

IF(N.GF.NSTOP) GO TO 20 
TEMT=TEMT+T 

IF(TEMT.GE ..16E-C8) GO TO 20 
NM 1 = N 
GC TO A 
20 KKK=KKK+ 1 

I I=KKK-KKK IN 
XP( I I 1 =TI ME 
YP ( I I 1 =-V ( N 1 
NM 1 = N 
N= 1 

TEMT=O.C 

IF (KKK .LT . 1493 1 GO TO 4 



100 



XP(1 >=0.0 
YP( 1 » = 10 . 

C PLOT THE DATA POINTS 

IF (I END. EQ. 1 ) GO TO 104 
IF( ICNT.GT .1 ) GO TO 1C2 

CALL DRAW ( 900, XP,YP t 1,0, LABEL 1, I T I TL E , 3. E- 07 , 3 . , 1 , 0 , 2 , 

10. 5. 5.1, L) 

WP IT El 6, 13 ) L 
13 FORMAT (//, IX, • L = ' , 110 ) 

GO TO 103 

102 I F ( I CNT . EQ . I END ) GO TO 101 

CALL DRAWl 900, XP,YP,2 ,0»LABEL1» IT I TLE , 3 . E-07 , 3 . , 1 , 0 , 2, 
10,5,5, 1 » L ) 

WR IT E ( 6 , 1 3 ) L 

103 I C NT = I CNT+ 1 
GO TO 100 

101 CALL DRAWl 900, XP,YP, 3,0, LAB ELI, ITITLE,3.F-07,3.» 1,0,2, 

10.5.5. 1 , L > 

WRITE16, 13) L 
GC TO 106 

104 CALL DRAWl 900 , XP , YP , 0 , 0 , LABEL1 , I T I TLE , 3. E-07 , 3 . , 1 , 0 , 2 , 
10,5,5, 1 , L ) 

WR ITE ( 6, 1 3 ) L 
196 RETURN 
END 



4. SOLUTION AND SENSITIVITY ANALYSIS PROGRAM FOR NONLINEAR 
DIODE EXAMPLE 



REAL*8 I T I TLE 

DIMENSION ITITLEl 12) .CURRl 600 > ,XP( 905 ) , V ( 600 > , YP ( 905 > , 
IX S If 600) ,S 1(905), XS21 600) ,S2( 905) ,XS3( 600) ,S3( 905) , 
2XS4160C) ,S4( 90 5) , SCURl 600 >,XS5(600),S5(905),XS6( 600) , 
3S6 ( 905 ) , TC UR ( 600 ) 

CALL CANCEL! 2) 

NST0P=3C0 . 

READ ( 5 , 12 ) LAB EL 1 
12 FORMAT ( A4 ) 

TI ME=G .0 

RFA0(5,3)T ,SI ,E,R,TR,TD,RADC,CD 

RE AD (5, 3) VI 1) , XP(1),YP(1), CURRl 1), SCURl 1) ,TCUR( 1) 

RE ADI 5,3) VM IN, SI MIN, S2 MI N , S3M I N , S4MI N , S 5 M IN , S6M I N 
RE AD (5, 3) XSll 1 ) . XS21 1 ) , XS31 1 ) ,XS4( 1) , XS51 1) ,XS6( 1 ) 

RF AD ( 5 , 3 ) Sill), 52(1), S3( 1),S4(1),S5( l),S6(l) 

3 FORMAT ( 8E1C. 2) 

TEMT=0 .0 

ALP= 1 • /. 02 6 
RSH=l.E+08 
G= 1 . /R+l ./ RSH 
UG=-E/R 
US10=C . 

US20=-E/R 
US 30=E /R 
US40=C . 

US 50=0 . 

US60=C . 

US 1N=C • 

US2N=-F/R 
US3N=E/P 
US 5N=C • 

US 1T=0 . 

US2T=US20+US2N 
US 3T=U S30+ US 3N 
US5T=0. 

GOLD=R ADC 
KKK= 590 
KKKI N= KKK- 1 
NM 1= 1 
N= 1 

4 IFITIME.LT.1.E-07) GO TO 1 
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GNFW=0 .0 
GO TO 2 

1 GNEW=RADC 

2 N=N+ 1 

T = 2 . *T 

IF(T.GE.l.E-lO) T= 1. E-10 
C0N1=S I*ALP*EXP( ALP*V< NM1 ) ) 

CGN2 = TD*C0N1 
CGN3=ALP*CCN2 

CGN4=S I*(EXP(ALP*V(NMD)-1.) 

DEN=CD+C0N2 

DEN2=DEN**2 

DEN3=DEN**3 ' 

F=- ( G* V( NM D + C0N4 ) /DEN 
B-l./DEN 

PF=- (G+COND/DEN + ( C0N3*( G*V(NMD +C0N4) ) /DEN2 
PB=-CO N3/D EN2 

F S 1=— ( (G+C0N1) *XS1(NMD )/DEN+( ( C0N4+G* V( NMD )*(C0N2 
1+CGN3*XSH NMD ) ) / D EN2 

F S2=- ( (G + CCND *XS2 (NM1 ) )/DEN+(CON3*XS2(NMl )*(G*V(NMD + 
1C0N4) ) /DEN 2 

F S3=- ( (G+C0N1 ) *XS 3 ( NM 1 ) -V ( NM D / R ) / DEN + ( CON 3*XS3(NMD* 
1<G*V(NMD+C0N4) ) /DEN2 

F S4= - ( (G+C0N1 ) *XS4(NM1 ) ) /DEN + ( CON 3*XS 4 ( NM 1 ) * ( G*V < NM1 ) + 
1C0N4) ) /0EN2 

FS5=-( (G+COND *XS 5 ( NM D ) /DEN+( { G*V ( NM1 ) + C0N4 I * ( CD + C0N3 
1*X S5 ( NM1 ) ) )/DEN2 

FS6=-( (G+CCND *XS6(NM1 ) )/DEN+(C0N3*XS6(NMl )* { G*V { NM1 ) + 
1C0N4 ) ) /DEN 2 

BS 1=- ( CON 2 +CGN 3*XS1(NMD J/DEN2 
BS2=-(C0N3*XS2(NMD )/DEN2 
BS 3=-( CON 3 *X S3 { NM 1 ) )/0EN2 
BS4=-(C0N3*XS4(NM1 ) )/DEN2 
B S 5=- ( CD+C GN3* XS 5 ( NM1 ) )/DEN2 
BS 6=- ( CON 3*X S6 ( NM 1 ) )/DEN2 

PFS1=- (G+C0N1 ) /DEN+(C0N3*(C0N4+G*V (NM1 ) ) )/DEN2 

P F S2=P FS 1 

PFS3 = P FS 1 

PF S4=PFS 1 

PFS5=PFS 1 

PF S6 = PFS 1 

P XF S 1=- ( AL P*C ON 1* XS1 ( NM1 ) ) /DEN + ( ( G+C0N1 ) * < C0N2+2 . *CON 3 
1*X S 1 ( NMD ) + (G*V(NMD+CON4)*(CON3 + ALP*CON3*XSl(NMD ) )/ 
2DEN2-( 2.*C0N3*(G*V(NM1 )+C0N4)*( CON2+CON3*X S 1 (NM1) ) )/ 
3DEN3 

PXFS2=-( ALP*C0N1*XS2( NM1) ) /DEN+ < ALP*C0N3*X S2 ( NM1 ) *<G* 
1V( NM1 ) +C0N4 ) +2 .*C0N3*XS2 ( NM1 ) *(G+C0N1 ) )/DEN2- ( 2.*C0N3* 
2C0N3*XS2(NM1 ) * ( G* V( NM D +C0N4 ) ) /DEN3 
PXFS3=-( ALP*C0N1*XS3( NMD-1. /R)/OEN + ( ALP*CON3*XS3( NM1 ) 
1*(G*V( NM1 ) +C0N4) +2.*C0N3*XS3(NM1 ) * ( G+CON 1 ) -C0N3* 
2V(NMD /R)/DEN2-( 2.*CON3*CON3*(G*V( NM1 ) +C0N4) *XS3( NM1 > ) 
3/DEN3 

PXFS4=-(ALP*C0N1*XS4(NMD ) /DEN+ ( AL P*C0N3 *X S4 ( NM1 ) *(G* 
IV (NM1) +C0N4)+2.*C0N3*XS4(NM1 )*( G+COND ) /DE N2- ( 2. *C 0N3* 
2C0N3*XS4(NM1 )* <G*V(NM1 >+C0N4) )/DEN3 
PXFS5=-( ALP*C0N1*XS5< NMD ) /DEN+ ( ( G + CON D *( CD + 2 .*C0N3* 

1 X S 5 ( NM 1 ) ) + { G*V ( NM 1 ) +C0N4 ) * ( ALP*C0N3*XS5( NMD ) )/DEN2-( 
22.*C0N3*(G*V(NM1 ) +C0N4 )*(CD+C0N3*XS5(NM1 ) ) )/DEN3 
PXFS6=-( ALP*C0N1*XS6( NMD ) /DEN+(ALP*C0N3*XS6(NM1 ) *(G* 

1 V ( NM1 ) +C0N4)+2 .*C0N3*XS6 (NMD* (G+CON 1 ) ) /DEN2- ( 2. *C0N3* 
2CCN3*XS6(NM1 )* (G*V(NM1 ) +C0N4) ) / DEN3 
PBS 1=-CGN3/DEN 2 
PBS2=PBS 1 
PBS3=P BS 1 
PBS4=P BS 1 
PBS5=PBS1 
PBS6=P BS 1 

PXBS l = -( C0N3*( 1 . + ALP *X S 1 ( NM 1 ) ) i /DEN 2+ ( 2. *C0N3* ( CON 2+ 
1C0N3*XSDNMD ) )/DEN3 

PXBS2=-< ALP*CON3*XS2( NMD ) /DEN2+ ( 2 . *CON3*CON3*XS2 ( NM1 ) 
D/DEN3 

PXBS3=-( ALP*CON3*XS3( NMD J /DEN2+ ( 2 . *C0N3*C0N3*XS3 ( NM 1) 
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1 )/DEN3 

PXBS4=-(ALP*C0N3#XS4(NM1) )/DEN2+(2.*C0N3*C0N3*XS4(NMl) 
1) /DEN3 

PXBS5=-(ALP*C0N3*XS5( NMD ) /DEN2+ ( 2 . *C0N3* ( CD+C0N3* 

1 XS 5 ( NM 1) ) ) /DEN3 

PXBS6=-( ALP*C0N3*XS6< NMD ) /DEN2M2.*CGN3*C0N3*XS6 I NM1 I 
1 )/DEN3 

9 CURR ( N ) = ( < 2.*TR-T )/(2 . *TR+T) ) *CURR ( NM 1 ) + ( T / ( 2 . *TR + T ) )* 
1 ( GOL D + GNEW ) 

SCUR (N )= ( ( 2 • *T R-T )/ (2. *TR + T ) ) *SCUR ( NM 1 ) + ( T /( 2.*TR+T) ) * 
1 ( GOLD+GNE W ) 

TCUR ( N )= ( ( 2 . * TR-T I + TCUR < NM 1 ) + T* ( CURR ( N ) +CURR ( NM1 ) 
1-GOLD-GNEW ) )/ ( 2.*TR+T ) 

UN=CURR(N)-E/R 
US4N=SCUR< NI 
US6N=T CUR ( N ) 

UT=UN+ UO 
US4T=US4N+US40 
US6T = USfcO+ US6N 

V< N) = V(NM1 ) + (T*( F+.5*B*UT D /( 1 .-T* ( PF+PB*UN) /2 . ) 
OELV=V(N)-V(NMl) 

VNCRM= ABS(DELV) 

IF( VNORM.GT.VMIN) GO TO 5 

XSD N) =XS1 (NM1 ) + ( T * ( FS 1 + . 5* ( DEL V* ( PXFS 1 + PB *US 1N+PXBS1* 
1UN)+B*US1T+BS1*UT) ) )/(l.-(T 1 MPFSl + PBSl*UN) )/2.) 

DF LS 1= XS 1 ( N)-XSl(NMl) 

SI NORM= SORT (DELV**2+DELS1* # 2) 

IF ( S1NGRM.GT.S1MI N ) GO TO 5 

XS2 ( N ) =XS2 ( NM 1 )+( T* ( FS 2 + . 5*( DEL V* ( PXF S 2+ PB *US2N+P XBS2* 
1UN) + B*US2T +BS 2 *UT ) ) ) / < I . - ( T* { PFS 2 +PBS 2 *UN ) )/2. ) 

DfcL S 2= XS2 ( N)-XS2( NMD 
S 2NO RM= SQRT( DELV**2+OELS2**2) 

IF (S2NORM.GT .S2MI N) GO TO 5 

XS3( N) = XS3(NM1 )+( T*(FS3+.5*( DEL V* ( P XF S 3+ PB *US3N+ PXBS3* 
1UN)+B*US3T+BS3*UT )D/(1.-(T*< PFS3+PBS3*UN) )/2. ) 

DE LS3= XS3( N) -X S3 ( NM1 ) 

S3N0RM= SORT ( DELV**2+DELS3**2 ) 

IFIS3NORM.GT.S3MIN) GO TO 5 

XS4( N)=XS4(NM1 ) + { T* ( F S4+ . 5* ( DE LV* ( PXFS4+PB*US4N+PXBS4* 
I UN )+B*US4T+BS4*UT D )/( l.-(T*(PFS4+PBS4*UN) )/2. ) 

DE LS4= XS4 ( N)-XS4( NMD 
S4N0RM= SORT (DELV*#2+DELS4**2) 

IFIS4N0RM.GT .S4MIN) GO TO 5 

XS5(N)=XS5(NMD+( T*(FS5+.5*< DELV*< PXFS5+PB*US5N+PXBS5* 
1UN ) + B*US5T +BS5*UT ) ) )/( l.-(T*(PFS5+PBS5*UN) )/2. I 
DELS5= XS5( N ) -X S5 ( NM1 ) 

S 5NORM= SQRT(DELV**2+DELS5**2) 

IF(S5NORM.GT .S5MIN) GO TO 5 

XS6( N) =XSMNM1 )+( T* ( F S6 + . 5* ( DE LV* < PXF S6 + PB*US6N + PX BS6* 
I UN ) + 6*US6T+BS6*UT) D /< l.-(T*< PF S6 + PB S6*UN ) )/ 2 . ) 
DELS6=XS6( N)-XS6 (NMD 
S6NORM= SORT (DELV#*2+DELS6#*2) 

IF (S6N0RM.GT.S6MIN) GO TO 5 
GO TO 7 
5 T=T/2. 

GO TO 9 
7 GOLD=GNEW 
TI ME=T IME + T 
UO=UN 
US 40=US4N 
US60=US6N 

IF (N.GE.NSTOP ) GO TO 20 
TEMT=T EMT + T 

IF(TEMT.GE ..15E-OB) GO TO 20 
NM 1=N 
GO TO 4 
20 KKK=KK K+ 1 

I I =KKK-KKK IN 
XP ( I I ) =TI ME 
YP ( I I ) =-V( N) 

SI ( I I )=XSi (N) 

S2( I I ) =XS2 (N) 
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on 



S 3 ( I I ) =XS3 ( N ) 

S4 ( I I >=XS4 (N) 

S5(II)=XS5(N) 

S6(I I )=XS6 IN ) 

WRITE (6, 14 ) TIME* CURR ( N ) ,V(N),XS1(N),XS2(N),XS3(N), 
1XS4(N) » XS5 ( N ) • X S6 ( N) » T »N 
14 FORMAT (1CE12.4* 15 ) 

NM1 = N 
N= 1 

TEMT=C .0 

IF (KKK.LT. 1493 ) GO TO 4 
PLOT THE DATA POINTS 
RE AD (5, 11) ITITLE 
11 FORMAT (6A8 ) 

CALL DRAW! 900 , XP , YP , 0 , 0 , LABEL 1 , 1 TI TLE , 3. E-07,3. , 

11.0. 2.2.5. 5. 1 , L ) 

WRITE (6, 13) L 

13 FORMAT (//, IX, * L= • , I 10) 

PLOT THE SENSITIVITY FUNCTIONS 
RE AD (5, 11) ITITLE 

CALL DR AW ( 900 , XP , S 1 , 0 , 0 , L AB EL 1 , I T1 TLE , 3. E-07 , 20. , 

10 .0. 2.2. 5. 5. 1, L) 

WRITE! 6, 13) L 

R E AD ( 5 , 1 1 ) ITITLE 

CALL DRAW ( 900, XP,S2, 0,0, LAB ELI, IT I TLE, 3. E-C7, 100 . , 

15.0. 2.2. 5. 5. 1, L) 

WP ITE( 6 , 13 ) L 

RE AD ( 5 , 1 1 ) ITITLE 

CALL DRAW (900, XP , S 3 , 0 , 0 , L AB EL 1 , I T I TLE , 3 . E- C7 , 1 00 . , 
10 ,0,2,2, 5, 5, 1, L) 

WRITE! 6, 13) L 
RF AD ( 5 , 1 1 ) ITITLE 

CALL DRAW! 900, XP, S4 , 0 , 0 , L AB EL 1 , I T I TL E , 3 . E- 07 , 1 00 . , 

10.0. 2.2.5, 5, 1 , L ) 

WRITE(6,13) L 

RE AD ( 5 , 1 1 ) ITITLE 

CALL DRAW 1900, XP, S 5 , C , C , L AB EL 1 , I T I TLE , 3. E-0 7 , 2 . , 

13.0. 2.2.5.5. 1 , L ) 



WR ITF ( 6, 1 3 ) L 
RE AD (5, 11 ) ITITLE 



CALL DRA W ( 900, XP, S6, 0,0, LABEL 1, I T I TL E , 3 . E- 07 , 1 0 . , 
11,0, 2,2, 5, 5, 1 , L ) 

WRITE (6, 13) L 
106 RETURN 
END 



5. TRANSISTOR BIASING PROGRAM 



CALL CANCEL! 2 ) 

C READ IN THE TRANSISTOR PARAMETERS 

RE AD (5, 2 ) RSHE ,CDE,TDE, SIE, ALPE, ALP I 
RE AD (5, 2 ) RSHC , C DC , TDC , S I C , ALPC , AL PN 

READ IN THE DESIRED BASE-EMITTER, BASE-COLLECTOR AND 
B I ASING VOLTAGES 
R E AD ( 5 » 2 ) E , X 1 , X2 
2 FORMAT (6E 1 C. 2 ) 

GC = (X1/PSHE-ALPI*SIC*(EXP( ALPC*X2)-1 . ) +S IE *( EXP! ALPE* 
1X1 )- 1. ) ) /( E+X2-X1 ) 

RC = 1 ./GC 

GP=(GC*(X1-E-X2 )-X2/RSHC+ALPN*SIE* ( EX P ( ALP E*X 1 ) - 1 • ) 
1-SIC*(EXP( ALPC*X2)-1. ) )/X2 
RB=1./GB 

C WRITE THE VALUES OF THE BIASING RESISTORS 

WR IT E ( 6, 5 ) RB , RC 
5 FORMAT! • 1' ,//2E15.5) 

RETURN 

END 
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6. SOLUTION PROGRAM FOR NONLINEAR TRANSISTOR EXAMPLE 



REAL*8 ITITLE 

DIMENSION F(2,1),DER(2,2),AI(2,2),TEMP(2,2), 

1VI9C5) ,XP(905) ,XS(2,1) , TE MP 1 I 2 , 2 ) , DER I NV I 2 , 2 ) , 

2BI 2, 2 ) ,U( 2 » 1 ) , BOER (2, 2) ,UT(2,1 ) , TEMP2 I 2 , 2 ) , UT EMP I 2 , 1 ) 
3, ITITLEI 12) 

CALL CANCELI2) 

NS=2 

NU=1 

READ! 5,20) ( ( A I ( I , J ) , J=1 ,NS ) , I =1 , NS ) 

20 FORMAT ( 8E 10. 2 ) 

PE AD ( 5 , 2 ) RSHE»CDE,RB , S I E , ALP E , AL P I 
R E AD ( 5 » 2 ) RSHC , CDC ,RC , S I C , ALPC , AL PN 
RE AD ( 5 , 2 ) E,TDC,TDE,T,TR,RADC 
2 FORMAT (6E10. 2) 

N= 1 
T I =0 . 

XM I N=. 01 
CM IN= 1 .F-Ofc 
TM AX= 1 • E— 0 8 
HCNT=1 . 

GE=1 ./RC+1 ,/RSHE 

GC=1./RC+1 ./RSHC+l./RB 

X 1 = . 6 

X2=-3. A 

VI 1 ) =X 1— X2 

XP 1 1 ) = T I 

XS(1,1)=X1 

XS (2 , 1 ) = X2 

UTEMPI 1, 1 ) =0. 

UTEMPI 2, 1 )=E 

6 IF ( T I .GE.l .5E-06 ) RADC=0. 

T= 2 . *T 

IF (T.GT.TMAX ) T = T M AX 
4 CCNE1=SIE*ALPE*EXP(ALPE*X1) 

C0NC1=SIC*ALPC*EXP(ALPC*X2 ) 

C0NF2=TDE*C0NE1 

CONC 2= TDC* CONC 1 

C0NE3=ALPE*C0NE2 

CCNC3= ALPC*CQNC2 

CGNE4=SIE*(EXP(ALPE*X1)-1.) 

CONC 4= SIC* (EXP(ALPC*X2)-1. ) 

DFNE=CDE+CCNE2 
DENC = CDC+CCNC2 

F ( 1, 1)=(-GE*X1+X2/RC+ALPI*C0NC4-C0NE4) /DENE 
F(2,l) =<-GC*X2+Xl/RC+ALPN*C0NE4-C0NC4)/DENC 
B ( 1, 11*1. /DENE 
B( 1,21=1./ (DENE*RC ) 

B ( 2 , 1 ) =C‘ . 

B(2,2)=-1./(DENC*RC) 

U ( 2 * 1 ) =E 

DERC 1, 1)=- (GE+C0NE1+C0NE3*F( 1,1) l/DENE 
DERI 1 , 2 ) = I l./RC+ALPI *C0NC 1 ) /DENE 
DERI 2, l) = t 1. /RC + ALPN*C0NE1 )/DENC 
DERI 2,2) =— (GC + CONC 1+C0NC3*F(2»1 ) )/DENC 
BDER(1,2)=C. 

BDER (2,1 )=C . 

8 U( 1, 1)=( I 2 • * TR— T ) *UTE MP ( 1 » 1 ) +2 . *T*RADC ) / I 2 .*TR+T ) 

BDER (1,1 ) = -(C0NE3*(U( 1 , 1 ) +E/RC ) ) / I DENE** 2 ) 

BDER (2 ,2 ) = (C0NC3*E )/ I RC* I DEMC**2 ) ) 

CALL ADD I BDER, DER , NS, NS, TEMPI ) 

Q=T/2 • 

CALL CONST I Q , TEMPI , NS , NS , TEMP ) 

CALL SUB I AI , TEMP, NS, NS, TEMPI) 

CALL RECIP I .000001, TEMPI, DERINV,KER, NS) 

Q=T 

CALL CONST I Q ,DER I NV , NS , NS , TE MP ) 

CALL ADD (U,UTEMP,NS, NU , UT ) 

Q= . 5 

CALL CONST I Q , UT , N S, NU , UT ) 
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noon 



CALL PROD (B,UT,NS,NS ,NU,TEMP2 ) 

CALL ADD ( TEMP2,F , NS , NU , TE MP1 ) 

CALL PROD (TEMP, TEMPI, NS, NS, NU.TEMP2) 

CALL ADD (TEMP2, XS,NS,NU,TEMP ) 

PROCEDURE GOOD FOR SMALL 
V E ( N ) -V E ( N- 1 ) 

VC< N)-VC( N-l ) 

DEFINE A NORM = SQRT(DELX1**2+DELX2**2) 

DELXI = TEMP( l , 1 > — X S ( 1,1) 

DELX2=TEMP (2,1 ) -X S ( 2 , 1 ) 
XN0RM=SQRT(0ELX1**2+DELX2**2 ) 

C NOR M= ABS ( U ( 1 , 1 ) -UT EMP ( 1 , 1 ) ) 

IF (CNURM.GT.CMIN) GO TO 5 
IF (XNORM.GT.XMIN) GO TO 5 
GO TO 7 
5 T=T/2. 

GO TO 8 
7 T I =T I + T 

X 1 =TEMP( 1,1) 

X2=T EM P ( 2, 1) 

0 = 1 . 

CALL CONST ( Q , U, N S ,NU , UTE MP ) 

CALL CONST ( Q, TEMP , NS , NU , XS ) 

TEMT=HCNT*l.E-08 

IF ( T I . LT .TEMT) GO TO 6 

HCNT=HCNT + 1 . 

N = N+ 1 
XP ( N ) = TI 
V ( N) =X 1-X2 

IF(N.LT.9C3) GO TO 6 
RE AD (5, 11 ) ITITLE 

11 FORMAT (6A8 ) 

READ(5»12) LAB EL 1 

12 FORMAT ( A4 ) 

CALL DRAW( 900, XP, V ,0 , 0 , L AB EL 1 , I T I TL E , 2 . E-06 , 1 . , 
10,0,2,2,5,5, 1 , L ) 

WRITE(6, 13) L 

13 FORMAT!//, IX,* L = ' ,110) 

106 RETURN 

END 
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